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The time and size distribution of the waves of topphngs in the Abehan sandpile model are expressed 
as the first arrival at the origin distribution for a scale invariant, time inhomogeneous Fokker-Planck 
equation. Assuming a linear conjecture for the the time inhomogeneity exponent as function of loop 
erased random walk (LERW) critical exponent, suggested by numerical results, this approach allows 
one to estimate the lower critical dimension of the model and the exact value of the critical exponent 
for LERW in three dimension. The avalanche size distribution in two dimensions is found to be the 
difference between two closed power laws. 
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The abelian sandpile model (ASM) was introduced by 
Bak, Tang, and Wiesenfeld as a minimal descrip- 
tion for natural phenomena characterized by intermittent 
time evolution through events called avalanches which 
have scale invariant properties. The model is defined on 
a hypercubic lattice whose sites can accommodate a vari- 
able, positive number of grains. With a uniform distribu- 
tion a site is chosen and its number of grains is increased 
by one. If its total number of grains exceeds a given crit- 
ical threshold Zmax, the nearest neighbor sites increase 
their number of grains by one and the initial site loses 
the corresponding number of grains. Then the newly up- 
dated sites are checked for stability until there are no 
more unstable states. This event is an avalanche and the 
analytical properties of its distribution is still an open 
question PJ3l0,@,o]. Analytical approaches rely on the 
algebraic properties of the toppling rules and the decom- 
position of avalanches into simpler events called waves 
which are related to spanning trees on a lattice 
A wave of topplings is simply obtained by restraining 
the initial site, from where an avalanche was initiated, to 
topple again only after all the other unstable sites have 
relaxed. Recently it was shown that the wave distribu- 
tion satisfies the finite size scaling ansatz, with critical 
exponents deduced from the equivalence between waves 
and spanning trees 

In this Letter we present numerical evidence that the 
critical exponents of the wave size and time distribution, 
which were deduced from geometrical considerations in 
[^,^1^ , can be related to the parameters of a scale invari- 
ant Fokker-Planck equation (FPE) in any dimension. In 
this way we make a connection between the geometrical 
properties of the waves and an evolution equation. As 
further confirmation of the validity of a FPE description 
for the abelian sandpile model, we show that the scaling 
behavior of the average number of unstable sites as func- 
tion of time is predicted by the FPE. Using the expres- 
sions for the critical exponents deduced in |^ together 
with those inferred form the FPE approach we are able 



to find the lower critical dimension of the abelian sand- 
pile model and, as an extra benefit, the exact value for 
the loop erased random walk (LERW) critical exponent 
V in three dimensions. In two dimensions, the case in 
which the distribution of the last wave is known |^], we 
can compute the asymptotic behavior of the avalanche 
distribution with the result that it has the form of a dif- 
ference between two close power laws. This has been 
previously proposed as an explanation for the failure of 
the finite size scaling approach 

The abelian sandpile model is, from its definition, a 
Markovian process whose states are specified by the lat- 
tice configuration. Once the initial point has been cho- 
sen randomly the dynamics of relaxation is deterministic, 
with the evolution determined by the initial configuration 
and the toppling rule. 

We consider a coarse grain description of the sand- 
pile evolution. Instead of a complete description, in- 
volving the number of grains at each site of the lat- 
tice, we use as variable the total number of unstable 
sites which exist at a given time, after an avalanche 
has started, irrespective of the configuration in which 
the sandpile is. This description is stochastic since the 
transition between two states with a given number of 
unstable sites depends on the configuration of the sand- 
pile which is now taken randomly. Let us consider for 
this process the evolution equation describing the tran- 
sition between states with different numbers of unsta- 
ble sites: Pit + \,n) = Y.„' W{t;n,n')P{t,n'); where 
W{n,n';t) is obtained by averaging over all the transi- 
tions between the configurations with the same number 
of unstable sites n, n' at a given time t. The transition 
probability W{t; n, n') may depend on time since the con- 
figuration of the unstable sites also depends on time; i.e. 
at the first step of a wave the unstable sites are some of 
the nearest neighbors of the initial site, eventually mov- 
ing away. In this formulation the wave is equivalent to a 
particle performing a discrete random walk on the posi- 
tive semi-axis with the transition probabilities depending 
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on the particle position. An wave event is a first arrival at 
the origin problem, for the wave stops when all the sites, 
except the initial one, are stable. This random walk is 
close to a diffusion process since the number of unsta- 
ble sites varies with bounded steps, x{t + 1) < 2Dx{t), 
where x{t) is the number of unstable sites and D is the 
lattice dimension. From this analogy we expect that the 
distribution of first arrival at the origin, that is the wave 
duration, has a power law distribution as it has for the 
first return distribution of the simple diffusion process. 

If we take the lattice size to infinity and the time unit 
to zero in an appropriate way the discrete Markov chain 
can be cast into a diffusion equation via the Kramer- 
Moyal expansion Here we shall not propose an ex- 
plicit way to construct the FPE for the ASM, instead 
we use the fact that in the stationary state the sand- 
pile is at criticality and we shall investigate the gen- 
eral FPE which yields critical behavior and has dif- 
fusion and drift coefhcients behaving similarly to the 
sandpile model. Generically the FPE has the form 
dtP= -d,j:[v{x,t)p] + {l/2)d'^^[D2ix,t)p], where p{x,t) is 
the probability density of the number of unstable sites, 
X is the number of unstable sites, t is the time since 
the wave has started. The drift coefficient v{x,t) and 
the diffusion coefficient D2{x,t) are obtained by tak- 
ing the continuum limit of the local first order moment 
~ jjt) and of the local second order mo- 

ment ^ 2;j)^iy(z, j; i) 10. Numerically we have 

found that the discrete diffusion coefficient £'2(2;, t) de- 
pends linearly on the number of unstable sites x. The 
slope depends on the dimension of the lattice but does 
not depend on the lattice size and the geometric condi- 
tion (bulk or boundary wave) ( see Fig. M . Also we have 
found that the discrete drift coefficient v(x) tends to a 
constant as the size of lattice grows. Fig. (|l|). The finite 
size effects affect the drift coefficient for the bulk waves 
at any value of the number of unstable sites, since the 
transition to states with larger number of unstable sites 
is smaller when the wave takes place near the boundary, 
while the statistics collects all the waves. 

The simplest Fokker-Planck equation satisfying the 
scale invariance assumption and the numerical behavior 
of the discrete coefficients D2(x,t) and v{x,t) is 

dtpix,t) = -d,[vt-y{x,t)] + dl,[D2xt-''p{x,t)] (1) 

where v,D2,a are constants. The initial condition for 
the above equation is p{x,t — to) = S{x — xq). Since 
we are interested in the time and size distribution for 
waves, which are first arrival events, we set an absorb- 
ing boundary condition at the origin p(x — 0, i) — 0, 
for the wave stops when the number of unstable sites, 
except for the initial one, is zero. The above differen- 
tial equation is invariant under a scale transformation 
X bx, t b^/^^~°'H. We observe that we can elimi- 
nate the parameter D2 by variable change x — > D2X and 
we change to v ^ v/ D2- 

Using a standard approach one can find the asymp- 



totic behavior of the first arrival at the origin forEq. (|): 
Pt(t)wt~^*, t>lwith 
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The second critical exponent we are interested in is the 
size distribution of waves. The size of the wave is the sum 
of the number of unstable sites until the wave stops; in 
the continuous formulation we have s{t) = j^^ dt'x{t') + 

We make the observation that the variable s is a mono- 
tonic function of time as s = x{t) > 0. The relation be- 
tween the two variables can be found using the average 
relation (s(<)) = {x{t)). Multiplying Eq. by x and 
integrating over x we obtain that {x{t)) « t ^ 1, 

where the average was normalized to the probability of 
surviving J dxp{x, t) until the moment t. Now we can use 
the variable s in the time distribution of waves. We have 
t ^^dt ~ s^-" ds^-" = s '^"ds for large t and s, where 
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This result can be checked easily by Monte Carlo simu- 
lation of a random walk constructed to be the discrete 
version of Eq. (|^) or equivalent ones. 

In Table | we show the values of the critical exponents 
Ta,Tt taken from using Eqs. we have computed 
the values of the parameters v and a. We observe that a 
has the same value for bulk and boundary waves and the 
time inhomogeneity disappears at the critical dimension, 
a = for Z? = 4. Thus the exponent a can be interpreted 
as a measure of the correlation among the unstable sites 
bellow critical dimension. 

Also, for D — 2 we note that a has the same numerical 
value as the critical exponent characterizing the decay of 
the autocorrelation function of waves found in One 
more hint that the exponent a is related to the corre- 
lation on the lattice is that it depends linearly on the 
erased random critical exponent v which has the values 

= 4/5 {D = 2),iy ^ 0.616 {D = 3), = 1/2 (D > 4) 
1^,^. Indeed one can see easily that the relation 
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holds exactly in D — 2 and D = A and has an error of 
0.003 for D — case in which the critical exponent v is 
only known from numerical simulation |ll|| . 

In the following we shall use Eq. H), which holds 
for both bulk and boundary waves, as a conjecture for 
further exploration of the Table I. The fact that a has 
the same value for both bulk and boundary waves can 
be checked numerically using the first two moments 
f^nii) = / x"p{x,t) of the solution of Eq. (|^). Inte- 
grating Eq. over x and using the absorbing boundary 
conditions at the origin, we obtain ■mi{t)/mo{t) « t^~°' 
which is independent of v. Numerical estimation of the 
above ratio in £' = 2,3, 4, presented in Fig. (g), shows an 
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excellent agreement with the predicted values (the error 
is less than 0.01). 

In 1^ it was shown that for waves the critical exponents 
for size and time distribution are given by 



Tt 1 + (d/ - (1 + cr))j^ 



(5) 



where is the fractal dimension of the wave, which has 
the values of the Euclidean dimension for D = 2,3,4 
and value 4 for higher dimension; = 1,0 for bulk and 
boundary wave respectively. 

Using the above proposed conjecture, Eq. (|^), and 
Eqs. (§§11), from (2 — a){Ta — 1) ~ n — 1 we obtain 
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This relation shows that that the minimum critical di- 
mension of the abelian sandpile model is 4/3 since the 
maximum value of z/ is 1. This result is in agreement 
with the observation that in _D = 1 the scaling in the 
abelian sandpile model breaks down |10|. 

An additional benefit of the above relation is that it 
gives the value of the critical exponent for the loop erased 
random walk in three dimension. Indeed if we put df — 
3, and keep in mind that = D for 13 < 4, we get 
1^0=3 = 8/13 = 0.615384... in perfect agreement with 
the numerical value found in pl| . In fact, assuming the 
conjecture Eq. (^ and identifying D = df, Eq. (^) yields 
a relation between the LERW critical exponent and space 
dimensionality. 

Ktitarev et al Q argue that all critical exponents of 
the abelian sandpile are determined by the critical expo- 
nent V. In order to complete this program we need the 
relation between the drift coefficient v and v and a rela- 
tion between the the drift coefficients v for the bulk and 
boundary waves. Using the Eq. (^) for Tt together with 
Eqs. (§,||) we obtain 
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We need one more relation to connect the coefficients 
a and v (bulk or boundary) with v. This can be obtained 
by using for example Eqs. (pP,H) from which we extract 
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V < 4/5 
> 4/5 
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We make the observation that the Eqs. (^,|7|,§D hold 
exactly in D — 3 if we use the deduced value v = 8/13, 
see Table I. At this point, we can conclude that we have 
found a self-consistent description of the critical proper- 
ties for the time and size distribution of waves in ASM. 
The exponent a captures the lattice correlation and the 
drift V controls the boundary condition of the wave, both 
of them being a function of the erased loop random walk 
critical exponent v through the Eqs. (§,0||). 



Now we are in the position to compute the asymptotic 
behavior for the avalanches in D = 2 using the FPE. 
In this description an avalanche is the sum of a random 
number of waves. The waves are statistically indepen- 
dent being recurrent events JTst . This assumption might 
appear to be in contradiction with analysis of jij but 
in fact in this approach the correlation is included al- 
ready in the time inhomogeneity. When a wave of size 
s touches the origin it has the the probability Pd{s) to 
die, thus also concluding the avalanche, or a new wave 
can start with the probability (1 —pd{s)). We choose the 
probability for the wave to die as Pdis) — s~^/^lns, so 
as to recover asymptotically the probability for the last 
wave: Pw{s)pd{s) — piw{s) ^ The avalanche 

distribution can then be written as Pa{s) — Pw{s)pd{s) -\- 
/o" ds'pw{s'){l -Pd{s'))pw{s - s')pd(s -s') + .... We can 
sum the previous series after a Laplace transform in s 
and we have 



Pa (A) 



Pito(A) 



1 - ((1 -pd)Pw)\ 



(9) 



Applying again the Tauberian theorem ||l4[| we find that 
asymptotically the avalanche size distribution behaves 
hke 



Pa{s) ~ Ci(slns) ^ + C2S - 



(10) 



This kind of behavior has been proposed previously in 
the literature (lO). The fact that l<ll/8<2 makes 
it difficult to obtain the 'pure' dominant behavior. From 
a numerical fit we obtain that C2/C1 ~ —0.25, there- 
fore Ci(slns)-i > C2S~^ for s > 10^. Thus, the 
FPE approach predicts that the avalanche distribution 
in the bulk must have the same asymptotic behavior as 
the waves for very large values of s, provided that the 
statistics excludes the avalanches which are affected by 
the boundary. 

In conclusion, we have used numerical hints to propose 
a FPE for the time and size distributions of the waves in 
the ASM. In this approach a wave is a first return event 
and the asymptotic properties of its distributions (time 
and size) are described by the first return probabilities 
of a time inhomogeneous FPE; the time inhomogeneity 
appears below the critical dimension D = 4. Further- 
more, this approach yields an analytical expression for 
the asymptotic behavior of the avalanche distribution in 
D = 2 which goes beyond the finite size scaling hypoth- 
esis and in agreement with recent results |^,|l^ . 

Using the relation for the critical exponents Tq, Tt de- 
duced in 1^ together with the the relations found through 
FPE approach and the conjecture (^ we the propose ex- 
plicit dependence of the critical exponents Tq, Tt of the 
critical exponent of LERW v (via a and v). A bonus 
of this approach is that it yields the value of the lower 
critical dimension, 4/3, for the ASM and the exact value 
of i/ 8/13 in D = 3 for LERW. 

The author thanks H. B. Geyer, F. Scholtz, L. Boon- 
aazier and A. van Biljon for useful discussions and 
H. B. Gcyer for a critical reading of the manuscript. 
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FIG. 1 

wave in sl)D = 2; L = 
waves and L = 1024 (□), 
b) 73 = 3; i = 128 (+) 
L = 128(D), L = 64(o) 



1024 (+), L = 512 (x) for bulk 
L = 512(o) for boundary waves; 
L = 64(x) for bulk waves and 
for boundary waves; c) D = 4; 
L = 32(+), L = 16 (x) for bulk waves and L = 32(D), 
L — 16(o) for boundary waves; d)the local second order mo- 
ment for bulk waves (+) {D = 2(L = 1024), D = 3{L = 128), 
D = 4(L = 32)) and for boundary waves D (D = 2(L = 512), 
D = 3{L = 64), D = 4(i = 16)). 
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a) D=2, L=512,1024 




b) D=3, L=64,128 
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FIG. 2. Finite size scaling for the moment ratio 
mio{t) = mi{t)/mo{t) in: D = 2(a), D = 3 (b), D = 4 
(c) with continuous hne for bulk waves and with dashed line 
for boundary waves. There is a good fit of the data with the 
exponents taken from Table (Q). The logarithmic correction 
to scaling at D = 4 has the same exponent as in Ref. 



TABLE I. The values for the critical exponent rt, Ta for 
the time and size distribution of waves taken from Q together 
with the values of the exponent a and the drift coefficient v 
computed from Eqs. (bUsI). The last line shows the values of 
the erased loop random walk critical exponent v. For D = 3 
we show in parenthesis the values computed with the exact 



value = 8 


/13. 
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